Analysis_01_prone_incidence_risk_reliability
Chad H Hochberg and Anna Barker
2025-08-22
#Specify Project Location and File Type
project_location <- '~/Library/CloudStorage/OneDrive-JohnsHopkins/Research Projects/CLIF/CLIF_Projects/ProneProjects/ProningEpi/ProneSevereARF_Output'
file_type <- 'csv'
#Create Sub Folders within Project Folder
# Check if the output directory exists; if not, create it
if (!dir.exists(paste0(project_location, "/summary_analysis"))) {
dir.create(paste0(project_location, "/summary_analysis"))
}df_agg <- read.csv(paste0(project_location, '/summary_output/aggregate_data.csv')) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
clif_hospital_description <- read.csv(paste0(project_location, '/clif_hospital_description.csv'))
#Merge Hospital Data Into DF AGG
df_agg <- df_agg |>
left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
#Define Hospital_Size Categorical Variable
mutate(icu_bed_cat=fcase(
Number_of_ICU_Beds<20, 0,
Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
Number_of_ICU_Beds>=50, 2
)) |>
## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
mutate(
icu_bed_cat = fcase(
Number_of_ICU_Beds < 20, 0L, # small
Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L, # medium
Number_of_ICU_Beds >= 50, 2L, # large
default = NA_integer_
)
) %>%
## 2. labelled factor ------------------------------------------------------
mutate(
icu_bed_cat = factor(
icu_bed_cat,
levels = 0:2, # make sure levels match
labels = c("small", "medium", "large")
)
)## Joining with `by = join_by(hospital_id)`#Do this Before Filtering
df_unique_hosp <- df_agg |> distinct(hospital_id, Hospital_Type)
cat('\nHow Many Academic v Community Hospitals are In THese Data')##
## How Many Academic v Community Hospitals are In THese Datatable(df_unique_hosp$Hospital_Type)##
## Academic Community
## 12 25df_hosptype_summary <- df_agg |>
group_by(Hospital_Type, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, Hospital_Type)## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.df_hosptype_summary## # A tibble: 6 × 6
## # Groups: Hospital_Type [2]
## Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
## <chr> <chr> <int> <int> <int> <dbl>
## 1 Academic COVID 12 1884 680 36.1
## 2 Community COVID 23 1085 643 59.3
## 3 Academic Post-COVID 11 1098 168 15.3
## 4 Community Post-COVID 23 498 150 30.1
## 5 Academic Pre-COVID 10 982 66 6.72
## 6 Community Pre-COVID 16 213 30 14.1#Hospital Size by Period
df_hospsize <- df_agg |>
group_by(study_period) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description')##
## Hospital Size Descriptiondf_hospsize## # A tibble: 3 × 4
## study_period small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 COVID 13 10 12
## 2 Post-COVID 12 11 11
## 3 Pre-COVID 8 8 10#Filter Out Hospitals with Less Than 10 Patients
df_agg <- df_agg |>
#Filter OUT HOspitals-Periods with < 10 observations
filter(n_patients>=10)
#Now After Filtering
cat('\nAfer Filtering Out Hospitals with >= 10 Patients \nHow Many Academic v Community Hospitals are In THese Data')##
## Afer Filtering Out Hospitals with >= 10 Patients
## How Many Academic v Community Hospitals are In THese Datatable(df_agg$Hospital_Type)##
## Academic Community
## 33 50df_hosptype_summary_post <- df_agg |>
group_by(Hospital_Type, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, Hospital_Type)## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.cat('\nProning Counts in Hospitals with >= 10 Patients')##
## Proning Counts in Hospitals with >= 10 Patientsdf_hosptype_summary## # A tibble: 6 × 6
## # Groups: Hospital_Type [2]
## Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
## <chr> <chr> <int> <int> <int> <dbl>
## 1 Academic COVID 12 1884 680 36.1
## 2 Community COVID 23 1085 643 59.3
## 3 Academic Post-COVID 11 1098 168 15.3
## 4 Community Post-COVID 23 498 150 30.1
## 5 Academic Pre-COVID 10 982 66 6.72
## 6 Community Pre-COVID 16 213 30 14.1#Hospital Size by Period
df_hospsize <- df_agg |>
group_by(study_period) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description')##
## Hospital Size Descriptiondf_hospsize## # A tibble: 3 × 4
## study_period small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 COVID 12 10 12
## 2 Post-COVID 7 11 11
## 3 Pre-COVID 3 7 10#Academic Size Summary
df_hospsize <- df_agg |>
group_by(Hospital_Type) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description by Hospital Type')##
## Hospital Size Description by Hospital Typedf_hospsize## # A tibble: 2 × 4
## Hospital_Type small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 Academic 0 3 30
## 2 Community 22 25 3#Proning By Hospital Size
df_hospsize <- df_agg |>
group_by(icu_bed_cat, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, icu_bed_cat)## `summarise()` has grouped output by 'icu_bed_cat'. You can override using the
## `.groups` argument.cat('\nProning Counts by Hospital Size and Period')##
## Proning Counts by Hospital Size and Perioddf_hospsize## # A tibble: 9 × 6
## # Groups: icu_bed_cat [3]
## icu_bed_cat study_period n_Hospitals n_patients n_proned percent_proned
## <fct> <chr> <int> <int> <int> <dbl>
## 1 small COVID 12 433 300 69.3
## 2 medium COVID 10 623 283 45.4
## 3 large COVID 12 1906 738 38.7
## 4 small Post-COVID 7 128 53 41.4
## 5 medium Post-COVID 11 335 77 23.0
## 6 large Post-COVID 11 1109 179 16.1
## 7 small Pre-COVID 3 43 6 14.0
## 8 medium Pre-COVID 7 156 23 14.7
## 9 large Pre-COVID 10 979 67 6.84#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned <- df_agg
df_agg_proned <- df_agg_unproned %>%
mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned <- df_agg_unproned |>
mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg <- bind_rows(df_agg_proned, df_agg_unproned)
#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
global_agg_byoutcome_period <- read.csv(paste0(project_location, '/summary_output/global_aggregate_byoutcome_period.csv')) |>
mutate(prone_log_odds = log(prone_rate_adjust / (1 - prone_rate_adjust))) |>
select(hospital_id, prone_12hour_outcome, prone_log_odds, study_period) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
#Merge Back with DF Agg
df_agg <- left_join(df_agg, global_agg_byoutcome_period) |>
#IF there were no patients proned can use the popensity for unproned patients by filling here
group_by(hospital_id, study_period) |>
arrange(hospital_id, study_period, prone_12hour_outcome) |>
fill(prone_log_odds, .direction = 'down') |>
ungroup()## Joining with `by = join_by(hospital_id, study_period, prone_12hour_outcome)`#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#Also Includes Hospital Type and Hospital Type and Period Interaction
#ICU Capacity is colinear with academic or community so have taken that out of the model
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% range
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Flat Priors for Hospital Type and Size
prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_full <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type +
(0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_full)
summary(agg_brm_period_full)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.91 0.12 0.69
## sd(study_periodPostMCOVID) 1.01 0.18 0.71
## sd(study_periodPreMCOVID) 1.00 0.26 0.58
## cor(study_periodCOVID,study_periodPostMCOVID) 0.77 0.11 0.51
## cor(study_periodCOVID,study_periodPreMCOVID) 0.60 0.21 0.12
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.52 0.23 0.01
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.18 1.00 1249
## sd(study_periodPostMCOVID) 1.40 1.00 1596
## sd(study_periodPreMCOVID) 1.59 1.00 1756
## cor(study_periodCOVID,study_periodPostMCOVID) 0.93 1.00 1683
## cor(study_periodCOVID,study_periodPreMCOVID) 0.93 1.00 1796
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.90 1.00 1623
## Tail_ESS
## sd(study_periodCOVID) 1858
## sd(study_periodPostMCOVID) 2226
## sd(study_periodPreMCOVID) 2729
## cor(study_periodCOVID,study_periodPostMCOVID) 2993
## cor(study_periodCOVID,study_periodPreMCOVID) 1940
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 1911
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## study_periodCOVID -0.11 0.25 -0.60
## study_periodPostMCOVID -1.20 0.30 -1.82
## study_periodPreMCOVID -2.10 0.33 -2.80
## Hospital_TypeCommunity 0.83 0.31 0.22
## study_periodPostMCOVID:Hospital_TypeCommunity 0.12 0.32 -0.49
## study_periodPreMCOVID:Hospital_TypeCommunity -0.16 0.49 -1.14
## u-95% CI Rhat Bulk_ESS Tail_ESS
## study_periodCOVID 0.40 1.01 656 1211
## study_periodPostMCOVID -0.63 1.00 942 1677
## study_periodPreMCOVID -1.50 1.00 1474 2067
## Hospital_TypeCommunity 1.43 1.00 650 1224
## study_periodPostMCOVID:Hospital_TypeCommunity 0.75 1.00 1721 2236
## study_periodPreMCOVID:Hospital_TypeCommunity 0.82 1.00 1998 1966
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## Elapsed Time: 2.163 seconds (Warm-up)
## Chain 2: 1.228 seconds (Sampling)
## Chain 2: 3.391 seconds (Total)
## Chain 2:
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.195 seconds (Warm-up)
## Chain 4: 2.235 seconds (Sampling)
## Chain 4: 4.43 seconds (Total)
## Chain 4:
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.402 seconds (Warm-up)
## Chain 1: 2.4 seconds (Sampling)
## Chain 1: 4.802 seconds (Total)
## Chain 1:#Now Generate Odds Ratios for Communtiy vs Academic Proning by Period
full_posterior <- as_draws_df(agg_brm_period_full) |>
mutate(pre_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPreMCOVID:Hospital_TypeCommunity`),
covid = exp(b_Hospital_TypeCommunity),
post_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPostMCOVID:Hospital_TypeCommunity`))
summary_posterior <- data.frame(
variable = c(
'Pre-COVID Community Hospital',
'COVID Community Hospital',
'Post-COVID Community Hospital'
),
median_or = c(median(full_posterior$pre_covid),
median(full_posterior$covid),
median(full_posterior$post_covid)),
conf.low = c(quantile(full_posterior$pre_covid, p=0.025),
quantile(full_posterior$covid, p=0.025),
quantile(full_posterior$post_covid, p=0.025)),
conf.high = c(quantile(full_posterior$pre_covid, p=0.975),
quantile(full_posterior$covid, p=0.975),
quantile(full_posterior$post_covid, p=0.975)))
print(summary_posterior)
write.csv(summary_posterior, paste0(project_location, '/summary_output/summary_posterior.csv'))####Predictions Command from Marginal Effects with 4000 posterior draws
reliability_adjust_full <- predictions(agg_brm_period_full,
type = 'response', ndraws=4000, re_formula = ~ (0 + study_period | hospital_id))
reliability_adjust_full <- as_tibble(reliability_adjust_full) |>
left_join(df_agg %>%
select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds', 'prone_log_odds', 'study_period', 'prone_12hour_outcome'),
by = c('hospital_id', 'study_period', 'prone_log_odds')) |>
#Need to Recombine Prone_12hour_outcome Estimates into 1 Per Hospital and Period
group_by(hospital_id, study_period) |>
mutate(
estimate_add=sum(estimate, na.rm=T),
ci_low_add=sum(conf.low, na.rm=T),
ci_hi_add=sum(conf.high, na.rm=T),
n_patients=sum(n_patients, na.rm=T)
) |>
arrange(hospital_id, study_period, -prone_12hour_outcome) |> #THis Arranges so Row Number 1 is the prone_12hour_outcome==1
filter(row_number()==1) |>
ungroup() |>
mutate(adjust_rate=estimate_add/n_patients) |>
mutate(ci_low=ci_low_add/n_patients) |>
mutate(ci_hi=ci_hi_add/n_patients) |>
mutate(period_rank=fcase(
study_period=='Pre-COVID', 0,
study_period=='COVID', 1,
study_period=='Post-COVID', 2
)) |>
#Orders for Graph
arrange(period_rank, adjust_rate) |>
mutate(hospital_rank=row_number()) |>
mutate(actual_rate=observed_prone/n_patients) |>
mutate(grouping='Actual Rate') |>
left_join(df_agg %>% select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds')) |>
distinct()## Warning in left_join(as_tibble(reliability_adjust_full), df_agg %>% select("hospital_id", : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 13 of `x` matches multiple rows in `y`.
## ℹ Row 13 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Joining with `by = join_by(hospital_id, Number_of_ICU_Beds)`## Warning in left_join(mutate(mutate(mutate(arrange(mutate(mutate(mutate(mutate(ungroup(filter(arrange(mutate(group_by(left_join(as_tibble(reliability_adjust_full), : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 151 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.break_points <- reliability_adjust_full |>
mutate(pre_covid_num=fcase(
row_number()!=1 & study_period=='Pre-COVID' & lead(study_period)=='COVID', hospital_rank),
post_covid_num=fcase(
study_period=='COVID' & lead(study_period)=='Post-COVID', hospital_rank
)) |>
summarise(
pre_covid_num=mean(pre_covid_num,na.rm=T)+0.5,
post_covid_num=mean(post_covid_num, na.rm=T)+0.5)
caterpillar_hospital_type <- ggplot(reliability_adjust_full,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.1),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals",
y = "Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
geom_vline(xintercept = break_points$pre_covid_num, linetype = 2, linewidth=0.3) +
annotate("text", x = 7, y = 0.95, label = 'Pre Pandemic', fontface = 'bold') +
geom_vline(xintercept = break_points$post_covid_num, linetype = 2, linewidth=0.3) +
annotate("text", x = 35.5, y = 0.95, label = 'Pandemic', fontface = 'bold') +
annotate("text", x = 68, y = 0.95, label = 'Post Pandemic', fontface = 'bold') +
theme_classic() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.5))
caterpillar_hospital_type
ggsave('caterpillar_hospital_byperiod_fully_adjusted.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
#With Actual Rate
#Add the Actual Value
with_rate <- caterpillar_hospital_type +
geom_point(aes(y=observed_prone/n_patients), color = 'maroon', alpha = 0.25)
with_rate
ggsave('caterpillar_fully_adjusted.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)#Stack Caterpillar Plots Vertically
vertical_data <- reliability_adjust_full |>
group_by(study_period) |>
arrange(study_period, adjust_rate) |>
mutate(hospital_rank=row_number()) |>
ungroup()
vertical <-ggplot(vertical_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.25),
strip.background = element_rect(fill = "grey90", color = NA),
strip.text = element_text(color = "black"),
legend.position = c(0.1, 0.90),
legend.background = element_rect(fill = "white", color = "grey"),
legend.title = element_text(size = 9)) +
facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')),
ncol=1,
scales='free_x')## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.ggsave('vertical_caterpillar.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
verticald <-ggplot(vertical_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.25),
strip.background = element_rect(fill = "grey90", color = NA),
strip.text = element_text(color = "black"),
legend.position = c(0.75, 0.85),
legend.background = element_rect(fill = "white", color = "grey"),
legend.title = element_text(size = 12),
legend.key.height = unit(1.5, "lines"),
legend.key.width = unit(2.2, 'lines')) +
facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')),
ncol=1)
ggsave('vertical_caterpillar_option_D.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
stacked_data <- vertical_data |>
group_by(study_period) |>
arrange(study_period, adjust_rate) |>
mutate(hospital_rank=if_else(
study_period=='COVID', row_number(), NA_real_
)) |>
group_by(hospital_id) |>
fill(hospital_rank, .direction = 'updown')
stacked <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, color = study_period)) +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= study_period), alpha = 0.4) +
geom_point(aes(color=study_period, shape=Hospital_Type), size = 2.5) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by COVID Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank())
ggsave('stacked_caterpillar.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).stacked2 <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, linetype=study_period), color = 'black') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) +
geom_point(aes(color=Hospital_Type, shape=study_period), size = 2) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Demuth") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by COVID Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
legend.position = c(0.15, 0.75),
legend.background = element_rect(fill = "white", color = "grey"))
ggsave('stacked_caterpillar2.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).stacked_data <- stacked_data |>
group_by(study_period) |>
mutate(period_n=sum(n_patients)) |>
ungroup()
stacked3 <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, alpha=period_n), color = 'black') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) + scale_alpha(guide="none")+
geom_point(aes(color=Hospital_Type, shape=study_period), size = 2.25) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Egypt") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Pandemic Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals",
shape="Study Period",
color= "Hospital Type"
) +
scale_shape_manual(values= c("circle", "triangle", "square"), labels=c("Pandemic", "Post-pandemic", "Pre-pandemic"))+
theme_classic() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
legend.position = c(0.15, 0.75),
legend.background = element_rect(fill = "white", color = "grey"))
ggsave('stacked_caterpillar3.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).posterior <- as_draws_df(agg_brm_period_full)
# Extract the standard deviations (SDs) for each study_period random effect
sd_pre <- posterior$`sd_hospital_id__study_periodPreMCOVID`
sd_covid <- posterior$`sd_hospital_id__study_periodCOVID`
sd_post <- posterior$`sd_hospital_id__study_periodPostMCOVID`
#3 Calculate MOR for Each Period
# MOR calculation function
calc_mor <- function(sd_samples) {
exp(sqrt(2) * sd_samples * qnorm(0.75)) # qnorm(0.75) ~ 0.674
}
# Compute MORs
mor_pre <- calc_mor(sd_pre)
mor_covid <- calc_mor(sd_covid)
mor_post <- calc_mor(sd_post)
# Differences Using The Latter Periods as Reference
diff_covid_pre <- mor_covid - mor_pre
diff_post_pre <- mor_post - mor_pre
diff_post_covid <- mor_post - mor_covid
# Posterior probabilities that the Difference in mOR is > 0 (Bayesian 'p-values')
prop_covid_gt_pre <- mean(diff_covid_pre > 0)
prop_post_gt_pre <- mean(diff_post_pre > 0)
prop_post_gt_covid <- mean(diff_post_covid > 0)
# Credible intervals
cr_covid_pre <- quantile(diff_covid_pre, c(0.025, 0.5, 0.975))
cr_post_pre <- quantile(diff_post_pre, c(0.025, 0.5, 0.975))
cr_post_covid <- quantile(diff_post_covid, c(0.025, 0.5, 0.975))
# Summarize
mor_summary <- data.frame(
Period = c("Pre-MCOVID", "COVID", "Post-MCOVID"),
MOR_Mean = c(mean(mor_pre), mean(mor_covid), mean(mor_post)),
MOR_Median = c(median(mor_pre), median(mor_covid), median(mor_post)),
MOR_2.5 = c(quantile(mor_pre, 0.025), quantile(mor_covid, 0.025), quantile(mor_post, 0.025)),
MOR_97.5 = c(quantile(mor_pre, 0.975), quantile(mor_covid, 0.975), quantile(mor_post, 0.975))
)
print(mor_summary)## Period MOR_Mean MOR_Median MOR_2.5 MOR_97.5
## 1 Pre-MCOVID 2.679054 2.514126 1.742355 4.564485
## 2 COVID 2.390461 2.355758 1.934036 3.075199
## 3 Post-MCOVID 2.655445 2.576878 1.968432 3.801425#Summarize the Differences in Median Odds Ratios by Period
diff_mor <- data.frame(
contrast = c('COVID-Pre', 'Post-COVID', 'Post-Pre'),
MOR_diff_mean = c(mean(diff_covid_pre), mean(diff_post_covid), mean(diff_post_pre)),
MOR_diff_median = c(median(diff_covid_pre), median(diff_post_covid), median(diff_post_pre)),
MOR_diff_2.5 = c(quantile(diff_covid_pre, 0.025), quantile(diff_post_covid, 0.025), quantile(diff_post_pre, 0.025)),
MOR_diff_97.5 = c(quantile(diff_covid_pre, 0.975), quantile(diff_post_covid, 0.975), quantile(diff_post_pre, 0.975)),
prob_diff_gt_0 = c(prop_covid_gt_pre, prop_post_gt_covid, prop_post_gt_pre)
)
print(diff_mor)## contrast MOR_diff_mean MOR_diff_median MOR_diff_2.5 MOR_diff_97.5
## 1 COVID-Pre -0.28859265 -0.1509414 -2.2549135 0.7904085
## 2 Post-COVID 0.26498433 0.2194309 -0.5086953 1.3141454
## 3 Post-Pre -0.02360833 0.0553715 -1.9437305 1.4517225
## prob_diff_gt_0
## 1 0.39250
## 2 0.71875
## 3 0.53575#Visualize
mor_df <- data.frame(
MOR = c(mor_covid, mor_pre, mor_post),
Period = factor(rep(c("COVID", "Pre-MCOVID", "Post-MCOVID"), each = length(mor_covid)))
)
ggplot(mor_df, aes(x = MOR, fill = Period)) +
geom_density(alpha = 0.5) +
labs(
title = "Period-Specific Median Odds Ratios (MOR \n (Fully Adjusted))",
x = "Median Odds Ratio (MOR)",
y = "Density"
) +
scale_x_continuous(seq(1,20, by=1), limits =c(1,20),
name='Distribution of Median Odds Ratos by Study Period') +
theme_minimal()
ggsave('mor_distribution_fully.pdf',
path = paste0(project_location, '/summary_output/graphs/'))## Saving 7 x 5 in image#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#Adjusted for COVID Status
#First Prepare Dataset
#Open Aggregate Data with COVID
df_agg_covid <- read.csv(paste0(project_location, '/summary_output/aggregate_data_wcovid.csv')) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
#Merge Hospital Data Into DF AGG
df_agg_covid <- df_agg_covid |>
left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
#Define Hospital_Size Categorical Variable
mutate(icu_bed_cat=fcase(
Number_of_ICU_Beds<20, 0,
Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
Number_of_ICU_Beds>=50, 2
)) |>
## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
mutate(
icu_bed_cat = fcase(
Number_of_ICU_Beds < 20, 0L, # small
Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L, # medium
Number_of_ICU_Beds >= 50, 2L, # large
default = NA_integer_
)
) %>%
## 2. labelled factor ------------------------------------------------------
mutate(
icu_bed_cat = factor(
icu_bed_cat,
levels = 0:2, # make sure levels match
labels = c("small", "medium", "large")
)
)## Joining with `by = join_by(hospital_id)`#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned_covid <- df_agg_covid
df_agg_proned_covid <- df_agg_unproned_covid %>%
mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned_covid <- df_agg_unproned_covid |>
mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg_covid <- bind_rows(df_agg_proned_covid, df_agg_unproned_covid) |>
#Only have the propensity at the period and outcome level so need to add period back here
mutate(study_period=fcase(
period_sarscov2=='COVID_SarsCov2Neg-Unk', 'COVID',
period_sarscov2=='COVID_SarsCov2Pos', 'COVID',
period_sarscov2=='Post-COVID_SarsCov2Neg-Unk', 'Post-COVID',
period_sarscov2=='Post-COVID_SarsCov2Pos', 'Post-COVID',
period_sarscov2=='Pre-COVID', 'Pre-COVID'
))
#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
#Merge Back with DF Agg COVID
df_agg_covid <- left_join(df_agg_covid, global_agg_byoutcome_period) |>
#IF there were no patients proned can use the popensity for unproned patients by filling here
group_by(hospital_id, period_sarscov2) |>
arrange(hospital_id, period_sarscov2, prone_12hour_outcome) |>
fill(prone_log_odds, .direction = 'down') |>
ungroup()## Joining with `by = join_by(hospital_id, prone_12hour_outcome, study_period)`set.seed(32284)
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% rnage
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PreMCOVID'),
#COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2Pos'),
prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2NegMUnk'),
#Post-COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2Pos'),
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2NegMUnk'),
#Flat Priors for Hospital Type and Size
prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_covid <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type +
(0 + period_sarscov2 | hospital_id)),
data = df_agg_covid,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start samplingplot(agg_brm_period_covid)
summary(agg_brm_period_covid)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type + (0 + period_sarscov2 | hospital_id)
## Data: df_agg_covid (Number of observations: 318)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 37)
## Estimate
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.94
## sd(period_sarscov2COVID_SarsCov2Pos) 0.91
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.96
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.60
## sd(period_sarscov2PreMCOVID) 0.96
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.86
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.79
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.69
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.55
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.52
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.61
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.48
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.56
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.43
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.37
## Est.Error
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.15
## sd(period_sarscov2COVID_SarsCov2Pos) 0.13
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.16
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.29
## sd(period_sarscov2PreMCOVID) 0.25
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.09
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.11
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.13
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.28
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.27
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.27
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.21
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.19
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.22
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.31
## l-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.69
## sd(period_sarscov2COVID_SarsCov2Pos) 0.69
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.69
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.06
## sd(period_sarscov2PreMCOVID) 0.57
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.64
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.51
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.38
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) -0.14
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) -0.13
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) -0.10
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.02
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.14
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) -0.03
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) -0.30
## u-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1.26
## sd(period_sarscov2COVID_SarsCov2Pos) 1.19
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.31
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1.24
## sd(period_sarscov2PreMCOVID) 1.52
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.97
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.95
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.89
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.92
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.91
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.94
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.84
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.88
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.81
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.86
## Rhat
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1.00
## sd(period_sarscov2COVID_SarsCov2Pos) 1.00
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## sd(period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1.00
## Bulk_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1740
## sd(period_sarscov2COVID_SarsCov2Pos) 2333
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2740
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 2148
## sd(period_sarscov2PreMCOVID) 2489
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 1576
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1916
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2309
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 4063
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 3853
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 3364
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 2899
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 3328
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 3592
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1615
## Tail_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 3077
## sd(period_sarscov2COVID_SarsCov2Pos) 2717
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2754
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1279
## sd(period_sarscov2PreMCOVID) 3011
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 2473
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2633
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2671
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 2710
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 2775
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 2474
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 3108
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 3586
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 3734
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 2482
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI
## period_sarscov2COVID_SarsCov2NegMUnk -0.98 0.23 -1.43 -0.53
## period_sarscov2COVID_SarsCov2Pos 0.47 0.23 0.02 0.90
## period_sarscov2PostMCOVID_SarsCov2NegMUnk -1.24 0.25 -1.73 -0.76
## period_sarscov2PostMCOVID_SarsCov2Pos -0.26 0.28 -0.81 0.26
## period_sarscov2PreMCOVID -2.28 0.30 -2.89 -1.73
## Hospital_TypeCommunity 0.83 0.27 0.31 1.33
## Rhat Bulk_ESS Tail_ESS
## period_sarscov2COVID_SarsCov2NegMUnk 1.00 991 1943
## period_sarscov2COVID_SarsCov2Pos 1.00 978 2073
## period_sarscov2PostMCOVID_SarsCov2NegMUnk 1.00 1303 2319
## period_sarscov2PostMCOVID_SarsCov2Pos 1.00 2277 3099
## period_sarscov2PreMCOVID 1.00 1810 2969
## Hospital_TypeCommunity 1.00 1023 1866
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).#Compare Differences in COVID vs Not COVID in COVID and POST COVID Period
post_covid <- as_draws_df(agg_brm_period_covid)
# Compute OR for the contrast
delta_covid_period <- post_covid$b_period_sarscov2COVID_SarsCov2Pos - post_covid$b_period_sarscov2COVID_SarsCov2NegMUnk
or <- exp(delta_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period \n', posterior_summary)## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period
## 4.241576 3.183033 5.564609# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2Pos - post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk
or <- exp(delta_post_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period \n', posterior_summary)##
## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period
## 2.655326 1.597988 4.562613#COVID Negative in COVID Period vs Pre-COVID
# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk - post_covid$b_period_sarscov2PreMCOVID
or <- exp(delta_post_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID \n', posterior_summary)##
## Odds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID
## 2.81339 1.60101 5.070301#ICU bed capacity Not Included as an additional covariable (in addition to hospital type becaues of co-linearity)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#ICU Capacity is colinear with academic or community so have taken that out of the model
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% range
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Flat Priors for Hospital Size
prior(normal(0,1), class = "b", coef = 'icu_bed_catlarge'),
prior(normal(0,1), class = "b", coef = 'icu_bed_catmedium'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_size <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + icu_bed_cat +
(0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_size)
summary(agg_brm_period_size)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + icu_bed_cat + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.95 0.13 0.72
## sd(study_periodPostMCOVID) 1.04 0.18 0.74
## sd(study_periodPreMCOVID) 1.01 0.26 0.59
## cor(study_periodCOVID,study_periodPostMCOVID) 0.80 0.10 0.56
## cor(study_periodCOVID,study_periodPreMCOVID) 0.64 0.19 0.20
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.57 0.22 0.09
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.23 1.00 1020
## sd(study_periodPostMCOVID) 1.44 1.00 1212
## sd(study_periodPreMCOVID) 1.60 1.00 1547
## cor(study_periodCOVID,study_periodPostMCOVID) 0.94 1.00 1564
## cor(study_periodCOVID,study_periodPreMCOVID) 0.93 1.00 1619
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.91 1.00 1549
## Tail_ESS
## sd(study_periodCOVID) 1687
## sd(study_periodPostMCOVID) 2334
## sd(study_periodPreMCOVID) 2129
## cor(study_periodCOVID,study_periodPostMCOVID) 2385
## cor(study_periodCOVID,study_periodPreMCOVID) 2136
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 2196
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID 0.67 0.25 0.15 1.13 1.01 804
## study_periodPostMCOVID -0.36 0.28 -0.94 0.15 1.01 842
## study_periodPreMCOVID -1.42 0.32 -2.08 -0.83 1.01 1099
## icu_bed_catmedium -0.30 0.36 -0.99 0.42 1.01 1016
## icu_bed_catlarge -0.53 0.33 -1.14 0.15 1.01 564
## Tail_ESS
## study_periodCOVID 1309
## study_periodPostMCOVID 951
## study_periodPreMCOVID 1780
## icu_bed_catmedium 1684
## icu_bed_catlarge 856
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## 00 [ 70%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.43 seconds (Warm-up)
## Chain 1: 1.399 seconds (Sampling)
## Chain 1: 3.829 seconds (Total)
## Chain 1:
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.515 seconds (Warm-up)
## Chain 2: 2.277 seconds (Sampling)
## Chain 2: 4.792 seconds (Total)
## Chain 2:
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.427 seconds (Warm-up)
## Chain 3: 2.567 seconds (Sampling)
## Chain 3: 4.994 seconds (Total)
## Chain 3:
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.524 seconds (Warm-up)
## Chain 4: 3.728 seconds (Sampling)
## Chain 4: 6.252 seconds (Total)
## Chain 4:#Now Generate Odds Ratios for Hospital Size Contrasts
full_posterior <- as_draws_df(agg_brm_period_size) |>
mutate(medium_v_small = exp(b_icu_bed_catmedium ),
large_v_small = exp(b_icu_bed_catlarge))
summary_posterior_size <- data.frame(
variable = c(
'Medium vs Small ICU Capacity',
'Large vs Small ICU Capacity'
),
median_or = c(median(full_posterior$medium_v_small),
median(full_posterior$large_v_small)),
conf.low = c(quantile(full_posterior$medium_v_small, p=0.025),
quantile(full_posterior$large_v_small, p=0.025)),
conf.high = c(quantile(full_posterior$medium_v_small, p=0.975),
quantile(full_posterior$large_v_small, p=0.975)))
print(summary_posterior_size)
write.csv(summary_posterior_size, paste0(project_location, '/summary_output/summary_posterior_size_analysis.csv'))#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
set.seed(32284)
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% rnage
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_strat <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_strat)
summary(agg_brm_period_strat)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.98 0.13 0.75
## sd(study_periodPostMCOVID) 1.09 0.17 0.78
## sd(study_periodPreMCOVID) 1.02 0.26 0.61
## cor(study_periodCOVID,study_periodPostMCOVID) 0.82 0.09 0.59
## cor(study_periodCOVID,study_periodPreMCOVID) 0.65 0.19 0.22
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.57 0.22 0.10
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.25 1.01 978
## sd(study_periodPostMCOVID) 1.46 1.01 1594
## sd(study_periodPreMCOVID) 1.61 1.00 1920
## cor(study_periodCOVID,study_periodPostMCOVID) 0.94 1.00 1968
## cor(study_periodCOVID,study_periodPreMCOVID) 0.94 1.00 1678
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.92 1.00 1063
## Tail_ESS
## sd(study_periodCOVID) 2105
## sd(study_periodPostMCOVID) 2546
## sd(study_periodPreMCOVID) 2812
## cor(study_periodCOVID,study_periodPostMCOVID) 2622
## cor(study_periodCOVID,study_periodPreMCOVID) 1814
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 1874
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID 0.39 0.17 0.05 0.72 1.00 660
## study_periodPostMCOVID -0.64 0.21 -1.06 -0.25 1.00 964
## study_periodPreMCOVID -1.74 0.26 -2.28 -1.27 1.00 1670
## Tail_ESS
## study_periodCOVID 874
## study_periodPostMCOVID 1536
## study_periodPreMCOVID 2356
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## ion: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.446 seconds (Warm-up)
## Chain 4: 1.113 seconds (Sampling)
## Chain 4: 2.559 seconds (Total)
## Chain 4:
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.489 seconds (Warm-up)
## Chain 2: 1.113 seconds (Sampling)
## Chain 2: 2.602 seconds (Total)
## Chain 2:
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.031 seconds (Warm-up)
## Chain 1: 1.196 seconds (Sampling)
## Chain 1: 3.227 seconds (Total)
## Chain 1:
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.703 seconds (Warm-up)
## Chain 3: 2.313 seconds (Sampling)
## Chain 3: 4.016 seconds (Total)
## Chain 3:#Posterior Draws
sensitivity_posterior <- as_draws_df(agg_brm_period_strat) |>
mutate(covid_v_pre_covid=exp(b_study_periodCOVID-b_study_periodPreMCOVID),
covid_v_post_covid=exp(b_study_periodCOVID-b_study_periodPostMCOVID),
pre_covid_v_covid=exp(b_study_periodPreMCOVID-b_study_periodCOVID),
post_covid_v_covid=exp(b_study_periodPostMCOVID-b_study_periodCOVID))
sensitivity_table <- data.frame(
variable = c('COVID_v_Pre-COVID',
'COVID_v_Post-COVID',
'Pre-COVID_v_COVID',
'Post-COVID_v_COVID'),
or = c(median(sensitivity_posterior$covid_v_pre_covid),
median(sensitivity_posterior$covid_v_post_covid),
median(sensitivity_posterior$pre_covid_v_covid),
median(sensitivity_posterior$post_covid_v_covid)),
conf.low=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.025),
quantile(sensitivity_posterior$covid_v_post_covid, p=0.025),
quantile(sensitivity_posterior$pre_covid_v_covid, p=0.025),
quantile(sensitivity_posterior$post_covid_v_covid, p=0.025)),
conf.high=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.975),
quantile(sensitivity_posterior$covid_v_post_covid, p=0.975),
quantile(sensitivity_posterior$pre_covid_v_covid, p=0.975),
quantile(sensitivity_posterior$post_covid_v_covid, p=0.975)))
print(sensitivity_table)
write.csv(sensitivity_table, paste0(project_location, '/summary_output/sensitivity_table.csv'))